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Abstract 

The phase structure of zero temperature twisted mass lattice QCD is investigated. 
We find strong metastabilities in the plaquette observable when the untwisted quark 
mass assumes positive or negative values. We provide interpretations of this phe- 
nomenon in terms of chiral symmetry breaking and the effective potential model of 
Sharpe and Singleton. 
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1 Introduction 

As a consequence of (soft) chiral symmetry breaking, nature has arranged itself such that 
three of the pseudo-scalar mesons are light, with masses around 140 MeV. This lightness 
of the pion mass becomes important also when we think of numerical simulations in lat- 
tice QCD. Approaching the "physical point", at which the pion mass assumes its value as 
measured in experiment, the algorithms used in lattice simulations suffer from a substantial 
slowing down ^ |2] which restricts present simulations to rather high and unphysical values 
of the quark mass. 

In addition to this slowing-down of the algorithms for Wilson fermions, the quark mass 
does not act as an infrared regulator allowing thus for the appearance of very small unphysical 
eigenvalues of the lattice Wilson-Dirac operator. These eigenvalues render the simulations 
more difficult and sometimes even impossible. 

Staggered fermions solve this problem but it is not clear how to use this approach to 
simulate Nf = 2 or odd number of flavours [3] . Overlap fermions |lj also solve the problem 
but they are computationally very demanding and, unless new algorithms are invented, they 
are very difficult to use for dynamical simulations. 

An elegant way out may be the use of so-called twisted mass fermions IE]- This 
formulation of lattice QCD is obtained when the Wilson term and the physical quark mass 
term are taken not parallel in flavour chiral space, but rotated by a relative twist angle u. 
If the Wilson term is given the usual form, such a chiral rotation leads to a twisted mass 
parameter /i, in addition to the standard Wilson quark mass ttlq ("untwisted" quark mass). 
Lattice QCD with a twisted mass was flrst employed for (9(a)-improved Wilson fermions 
with the nice feature that the improvement coefficients and the renormalisation constants 
are the same as for (9(a)-improved Wilson fermions without twisted mass and hence they 
did not need to be recalculated [Zj. The main advantage of the twisted mass fermions is 
that the twisted quark mass provides a natural infrared cut-off and avoids problems with 
accidental small eigenvalues rendering therefore the simulations safe. Of course, the slowing 
down of the algorithms when approaching small quark masses will remain, although it is 
expected to be less severe. 

Later on it was realized that a full (9(a)-improvement of correlation functions can be 
obtained by using the twisted mass alone without additional improvement terms when, as a 
special case, tuq is set to the critical value merit and the above mentioned twist angle is equal 
to to" = 7r/2 [H]. In this way the demanding computation of many improvement coefficients 
can be avoided rendering the simulations much easier both from a conceptual as well as from 
a practical point of view. 

The Wilson twisted mass formulation has been tested numerically in the quenched ap- 



proximation already jHj. The results are very encouraging. The 0{a) corrections appear 
indeed to be cancelled and even higher order effects seem to be small, at least for the quan- 
tities and the value of the quark mass considered in ref . [5] . 

A word of caution has to be added at this point. Although, as mentioned above, the 
twisted quark mass can be decreased towards zero without simulations breaking down due 
to exceptional configurations, there is an important interplay between the lattice cut-off, 
A = a~^ with a the lattice spacing, and the quark mass rriq (see equation (fTUI) below). In the 
continuum in the presence of spontaneous chiral symmetry breaking the chiral symmetry is 
not realized a la Wigner and, as the quark mass goes to zero, the chiral phase of the vacuum 
is driven by the phase of the quark mass term. The same must be true on the lattice, thus 
the scaling limit a — *> should be taken before letting nig ^ 0. As a result, taking the chiral 
limit is a numerically delicate matter. 

In order to ensure in practice that on the lattice the chiral phase of the vacuum is 
determined by the quark mass term, proportional to rriq, and not by the Wilson term, the 
lattice parameters should satisfy the order of magnitude inequality [Hj 

"^g^QCD > aAqcD • (1) 

This same condition emerges from many different corners of the lattice theory when the 
physical world is approached. A very simple argument leading to the bound (0) is obtained 
by comparing the magnitude of the critical Wilson term to that of the quark mass term and 
requiring the first to be negligibly small compared to the second one, in order to be sure 
that lattice physics matches the requirements of the continuum theory. From the order of 
magnitude inequality a(AQCD)^ *^ "^^(Aqcd)^, one immediately gets the condition (jT)). It 
is important to observe, however, that the less restrictive condition 
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Aqcd > («Aqcd)' (2) 



may be sufficient if one is dealing with (9(a)-improved quantities. 

It should be remarked that, since aAqco can be (non-perturbatively) expressed in terms 
of Qq, eqs. dH) and Q are actually (order of magnitude) conditions for the values of the 
dimensionless bare lattice parameters arriq and Qq. Contact with dimensionful quantities can 
be made by comparing simulation data with physical inputs. 

What is in practice important is to know for which range of the bare lattice parameters 
one can avoid troubles from chiral breaking cutoff effects, even if parametrically of order 
a^ or higher. This issue has to be settled by numerical investigations aimed at establishing 
both the structure of the phase diagram of the lattice model in study and the size of residual 
scaling violations on the physical observables. 



In this perspective, twisted mass fermions offer a unique opportunity to explore the phase 
diagram of Wilson fermions. By fixing the twisted mass parameter /i, one may vary {mQ — 
TTicrit) from positive to negative values. In this way, the phase diagram of zero temperature 
lattice QCD can be explored. It should be emphasized that, on large lattices, such an 
investigation would be very difficult without having yU 7^ 0, since else the algorithms would 
slow down dramatically approaching the critical quark mass. 

In this work we have performed simulations to explore the phase diagram of zero tem- 
perature QCD. As it will be shown in the following, we find strong metastabilities in the 
plaquette expectation value. We determined in both metastable branches a number of quan- 
tities such as the (untwisted) PCAC quark mass and pseudo-scalar meson masses. The 
results presented in this paper are obtained at only one value oi (3 = 5.2, with (3 related 
to the bare gauge coupling g^ hj P = G/g^. Since the value oi P = 5.2 corresponds to a 
rather coarse value of the lattice spacing (a ~ 0.16 fm) our work can only be considered 
as a starting point for a more detailed investigation of the phase diagram of lattice QCD. 
In particular, the /5-dependence of the strength of the observed metastabilities has to be 
determined. We believe that a qualitative and even quantitative understanding of the phase 
diagram is a necessary prerequisite for phenomenologically relevant numerical simulations. 

The paper is organized as follows. In section |21 we introduce Wilson twisted mass fermions 
and give our notations. This is followed by a short discussion of the algorithms used. In 
section El we provide our evidence for metastabilities by hysteresis effects and long living 
metastable states. There, we also show results for a selected set of physical quantities. 
In section IH we give a possible interpretation of these results in terms of chiral symmetry 
breaking and the Sharpe-Singleton effective potential model. We conclude finally in sectional 
In the appendix some details of the applied update algorithms are explained. 

2 Lattice action and basic variables 

2.1 Lattice action 

Let us start by writing the Wilson tmQCD action as 

S[U,x,x] = xiD[U]+mo + iJ,i-fr^Ts)x- (3) 



In eq. (jHl) rrio is the quark mass parameter and fj, is the twisted quark mass parameter. The 
operator D[U] is given by 

xD[U]x = a'J2\-x{x)x{x) (4) 
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- ^X(a;)^ {U{x, fi){r + -f^)x{x + afi) + U^x - afi, ^){r - -ff,)x{x - afi) ^ 



with r the Wilson parameter which will be set to r = 1 in our simulations. 

The action as it stands in eq. (jHl) can, of course, be studied in the full parameter space 
{mo, /i). A special case arises, however, when niQ is tuned towards a critical bare quark mass 
rrtcrit- In such, and only in such a situation all physical quantities are, or can easily be, 0{a) 
improved. It is hence natural to rewrite 

mo = merit + fh (5) 

with fh an offset quark mass. The values of merit need only to be known with 0{a) accuracy 
[Hj and can be, for instance, taken from the pure Wilson theory at /i = 0. 

For standard Wilson fermions usually the hopping parameter representation is taken in 
the numerical simulations. This representation is easily obtained from eq. Q by a rescaling 
of the fields 

a^/^ ' a^/'^ ' 2amo + 8r 

We then obtain the form of the action that is actually used in our simulations 

S[x.X.U] = ^|x(a;)(^l + 2m/ifi:75r3Jx(a;) (7) 

X 

4 

- i^x{x)^ {U{x, n){r + -f^)x{x + aft) + U^ {x - aft, fi){r - -f^)x{x - aft) ^ 

Although not needed for the discussion of the numerical data presented below, we give 
for completeness here the action in the so-called physical basis. This action is obtained by 
introducing new fields ip{x) and ^(x) which are related to the fields in eq. Q by a chiral 
transformation 

ip{x) = e*^'^5^3^(x) = (^cos - + i-f5T3 sin - j x{x) , 

i>{x) = x(a;)e*^'^'''' = X(a^) (cos - + 275X3 sin- j . (8) 



The action then reads 



S[ip,^p,U] = a^^ jmg^(x)V^(x)-— ^(x)e-^'^^^"^x 

X ^ 
4 

y irU{x, fi)ip{x + afi) + rU'{x — ap,, fi)ilj{x — afi) j — (2amci.it + 9>r)%l){x) 



1 -, 



-^^(^) X^ (yi.^^ ^^)l^li^{x + afi) -U\x- afi, ^i)i^,i){x - a/i) j > (9) 



where we have identified 

niq cos a; = uiq — merit = ""^5 ''^g sina; = /i . (10) 

2.2 Simulation algorithms 

In our numerical simulations we used two different optimized updating algorithms for pro- 
ducing samples of gauge configurations: the Hybrid Monte Carlo (HMC) algorithm with up 
to three pseudo-fermion fields as suggested in fOlE] and the two-step multi-boson (TSMB) 
algorithm ^2]- 

In the standard HMC algorithm we used even-odd preconditioning, which in presence of 
a twisted mass is only a slight modification of the standard preconditioning technique |13j . 
We give the relevant equations in appendix lA.il of this paper. As a subsequent improvement 
of the algorithm, we implemented the idea of ref. fH] and used shifted fermion matrices to 
"precondition" the original fermion matrix. These shifted matrices are treated by introducing 
additional pseudo-fermion fields. In the shifted fermion matrix we simply used larger values 
of the twisted mass parameter than the value of /i that is to be simulated. Using two 
pseudo-fermion fields we experienced a substantial improvement of the HMC algorithm by 
at least a factor of two. The addition of a third pseudo-fermion field gave only another 
10-20% improvement. Again we list the relevant equations, how the shifted matrices are 
implemented, in appendix lA.il As a further algorithmic trick we used the Sexton- Weingarten 
leap-frog integrator as proposed in ref. [T^ . 

Our alternative algorithm, the TSMB algorithm ^2], is based on the multi-boson repre- 
sentation of the fermion determinant |T3] . Optimized polynomial approximations are used, 
both in the first update step and in the second global accept-reject correction step, for repro- 
ducing the dynamical effect of fermions on the gauge field. We apply high order least-square 
optimization and obtain the necessary polynomials using high precision arithmetics |16j . 
Concerning the optimization of TSMB for QCD see, for instance, |17j . 

A useful improvement of the TSMB update algorithms can again be achieved by even-odd 
preconditioning. This can be implemented in TSMB for twisted mass quarks along the lines 



of ref. PHj- For the even-odd preconditioning of the TSMB update the flavour indices of the 
quark fields have to be kept. This means that the multi-boson fields have 24 components 
per lattice site (2 for flavour, 3 for colour and 4 for Dirac spinor indices). Correspondingly, 
the polynomials are approximating the function x~2 as in the case of a single Dirac flavour 
with untwisted quark mass. We give some more details of our even-odd implementation of 
the TSMB algorithm in appendix lA. 21 

In the region of light quarks an important part of the numerical effort has to be spent 
on equilibrating the gauge configuration in a new simulation point. This is particularly 
relevant in studies of the phase structure where many different points in the parameter 
space have to be investigated. In case of TSMB the equilibration time is substantially 
longer than the autocorrelation of relevant physical quantities in equilibrium: on our lattices 
equilibration can take ten or more times the autocorrelation time of the plaquette observable. 
The autocorrelation times in equilibrium themselves are similar but most of the time by 
factors of 2-3 shorter with our twisted quark masses than with untwisted quark masses of 
similar magnitude. For an approximate formula of the computational cost see ref. [TH] . 

The use of two different optimized update algorithms was very helpful in checking our 
results. We did not try to obtain a precise performance comparison. Qualitatively, we did not 
see a noticeable difference in the speed once equilibrium was reached, but the HMC algorithm 
with multiple pseudofermion fields (MPHMC) turned out to be faster in the equilibration 
process. In particular, crossing the transition region below and above the metastability 
region is faster with MPHMC. Nevertheless, the extension in n of the metastability region 
is the same with both algorithms. 

The data used for preparing the figures in this publication were obtained with MPHMC, 
except for the upper four panels in fig. |2l which were obtained with TSMB. The thermal 
cycles in fig. ^ were only run with MPHMC. In the other figures the results of the TSMB 
runs, whenever performed, were always consistent within errors with the shown MPHMC 
results. 

3 Numerical results 

In this section we give our numerical evidence for the phenomenon of metastability mentioned 
in the introduction. As a first step and for an orientation we have investigated thermal cycles 
in the hopping parameter n. We then discuss metastable states in the plaquette expectation 
value. Finally we determine quantities such as the pion mass and the untwisted PCAC quark 
mass in the metastable branches in order to obtain a picture of the physical properties in 
the different states. In most cases we perform the simulations at a twisted mass a/x = 0.01 



but in a few cases we also put a/x = 0, which is possible on the lattice sizes we consider. 

3.1 Thermal cycles 

We started our investigation of the phase diagram of zero temperature lattice QCD by 
performing thermal cycles in k, while keeping fixed /3 = 5.2 and the value of the twisted mass 
parameter afi. These cycles are performed such that a starting value of Kgtart is chosen and 
then K is incremented, without performing further intermediate thermalization sweeps, until 
a final value of Kfinai is reached. At this point the procedure is reversed and k is decremented 
until the starting value Kgtart is obtained back. At each value of k 150 configurations are 
produced and averaged over. 

In fig. ^ we show three such thermal cycles, performed at afi = 0, a/i = 0.01 and afi = 0.1 
from bottom to top. In the cycles signs of hysteresis effects can be seen for a/x = and 
afi = 0.01 while for the largest value of a/i = 0.1 such effects are hardly visible. Hysteresis 
effects in thermal cycles may be signs of the existence of a first order phase transition. 
However, they should only be taken as first indications. Nevertheless, they provide most 
useful hints for further studies to search for metastable states. 



3.2 Metastability 

Guided by the results from the thermal cycles, we next performed simulations at fixed values 
of afi and k, starting with ordered and disordered configurations, staying again a.t (3 = 5.2. 
In fig. |21we show the Monte Carlo time evolution of the plaquette expectation value, in most 
cases on a 12'^ x 24 lattice. For several values of k we find coexisting branches with different 
average values of the plaquette. The gap (the "latent heat") appears to be rather large. At 
K = 0.1717 we show the history of the plaquette expectation value also on a larger (16'^ x 32) 
lattice. It seems that the gap in the plaquette expectation value does not depend much 
on the lattice size, suggesting that the metastability we observe here is not a finite volume 
effect. In most cases the twisted mass is afi = 0.01, except for the picture left in the bottom 
line where it is a/x = 0. 

The lifetime of a metastable state, i.e. the time before a tunneling to the stable branch 
occurs, depends on the algorithm used. In fact, one may wonder, whether the appearance 
of the metastable states seen in fig. |21 may not be purely an artefact of our algorithms. We 
cannot completely exclude this possibility but we believe it is very unlikely: we employed two 
very different kinds of algorithms in our simulations as explained in section 12.21 We observe 
the metastable states with both of them. We also inter-changed configurations between the 
two algorithms: a configuration generated with the algorithm A was iterated further with 
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Figure 1: Thermal cycles in k, on 8^ x 16 lattices at /5 = 5.2. The 
plaquette expectation value is shown for: afi = 0.1 (A); afi = 0.01 
(B); a/i = (C). The triangles (V) refer to increasing /t-values, the 
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algorithm B and vice versa. We find that in such situations the plaquette expectation value 
remains in the state where it has been before the interchange of configurations took place. 
In addition, as we shall see below, the two states can be characterized by well defined and 
markedly different values of basic quantities. We therefore conclude that the metastable 
states are a generic phenomenon of lattice QCD in our formulation. 

3.3 Pion and quark masses 

By selecting separately configurations with high and with low plaquette expectation value, 
we measured the pion mass and the untwisted PCAC quark mass to study the physical 
properties in the two metastable states. 

We obtained the pseudo-scalar ( "pion" ) mass from suitable correlation functions. These 
are constructed from the standard composite fields defined in terms of the fields ip and ijj in 
eq. dHJ: 



S\x) - 


= '4){x)ip{x), 


P"{x) -- 


= H^hbY^ix) 


^^(^) = 


= '^{xh^i^Y'^ix), 


V^{x) 


= i^ix)'jf,^^{x 



Here Tq, a = 1, 2, 3 are the usual Pauli- matrices in isospin space. The corresponding compos- 
ite fields in terms of the quark fields x and x of eq. (jH)) are then given by the transformation 
in eq. (jH)). For instance, for a = 1,2 ("charged pions") the pseudo-scalar density has the 
same form in the x-basis as in the ■i/'-basis. Therefore, the mass of the charged pions can 
be extracted from correlators in the x-basis in the usual way. The charged axial vector and 
vector currents are rotated into each other by the angle uj in such a way that at cj = n/2 
they are interchanged. (For more details see the literature, e.g. |S1E]-) 

Besides the pion mass, we measured the PCAC quark mass from the axial vector current 
in the x-basis: 



PCAC ^ {9*^X1 ,.75^ x{x) 0^{y)) 
'"^ ~ 2(x^75X(x) OT(y)) 



PCAC ^ n::m:^^^^^^^av;^v^:^^^ _ ^^2) 



Here O^ is a suitable operator that we have chosen to be the pseudo-scalar density O^ = 
X^l5X{x) , 9* is the lattice backward derivative defined as usual and r"^ = ri ±2x2. One can 
show that in the limit a ^ the quantity mF'^^^ is asymptotically proportional, through 
finite renormalization constants, to rh. 

In fig.iniwe show the pion mass squared in lattice units as function of (2k)~^. We observe 
that the pion mass is rather large and the most striking effect in the graph is that it can 
have two different values at the same k. If we consider the quark mass mF.^'^^ in fig. |3J we 
see that in the states with a low plaquette expectation value the mass is positive while for 
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high values of the plaquette expectation it is negative. These quark masses with opposite 
sign coexist for some values of k. Plotting the pion mass versus m^'^^^ one obtains fig. El 
Figs. 12111] clearly reveal that for small enough fi metastabilities show up in the quantities 



we have investigated, such as m^, m\,^ ^ and the average plaquette, if uiq is close to its 
critical value. What "small enough /i" means is likely to change with (3. Simulations at 
larger values of (3 are in progress. As a matter of fact, when rriQ is significantly larger 
(smaller) than merit we find m^^^^ to be positive (negative) and no signal of metastabilities. 
The remark that metastabilities take place for rriQ close to its critical value will be important 
both in sect. 14.11 to understand why they affect also a purely gluonic observable such as 
the plaquette and in sect. 14. 3| where it leads to a plausible explanation of the observed 
metastability phenomena in terms of spectral properties of the lattice tmQCD Dirac matrix 
(suppression of the "eigenvalue cloud crossing" phenomenon by the fermionic determinant). 
The remarks in sect. 14.11 may provide further insight also on the similar metastability 
phenomena reported in [SDj for the Nf = 3 untwisted Wilson theory and on the reason 
why they "disappear" when changing the gluonic action or details, e.g. csw-'v^lne, of the 
fermionic action. A possible reason is that these changes might shift the range of niQ where 
metastabilities appear to values where no data are yet available. 



4 Physical interpretation 

The observed strong metastabilities discussed in the previous section clearly suggest that we 
are working either directly at a first order phase transition or at least very close to it such 
that we see the remnants of a close-by first order phase transition. With the present data we 
cannot really differentiate between these two scenarios and in the following we will therefore 
discuss both of them. 

4.1 Jump in the plaquette and chiral symmetry breaking 

Generally speaking, a jump in the plaquette as seen in our data can arise owing to the lack of 
chiral symmetry for chirally non-invariant formulations of lattice QCD. The argument relies 
on the key observation that, when working with chirally twisted Wilson fermions, there are 
two distinct sources of chirality breaking. The first source is the combination of the untwisted 
Wilson and mass terms 



xM[U]x = a 



^ixix)y— + mojxix) 



r 
2a 



^(^) X] (^*^^' '")^^^ + ^'"^ -U\x- afi, fi)x{x - afi)j > . (13) 
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Figure 3: The pion mass squared in lattice units on two lattice 
sizes measured separately on configurations in the two metastable 
states. These runs were made a.t /3 = 5.2 and afi = 0.01. 
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The second source of chirality breaking is the twisted mass term fiX'^'l^'^sX- As pointed out 
in section 121 one may trade the bare parameters rrio and /i in eq. (jH)) for the equivalent bare 
parameters nig and u of eq. ©• The latter are best suited to discuss the connection with 
continuum QCD physics, as u is an unphysical parameter, while rrig represents the bare quark 
mass. Assuming spontaneous chiral symmetry breaking in infinite volume, the pion mass 
squared is expected to vanish linearly in nig (up to lattice artifacts) as nig -^ 0. Moreover in 
the continuum limit the physical scalar condensate is expected to show a discontinuity and 
changes sign as rrig passes through zero: 

lim ([^V^]r) = - lim ([^^]r,)^0, (14) 

rrtq— >0+ mq— >0~ 

where by [^V^Jr we mean the appropriately subtracted and renormalized scalar density. We 
recall that f or uj ^ this is a non-trivial linear combination of XX^ X'^l5'^3X ^^^ the constant 
field (see below for details). 

In order to make contact with the observed metastability phenomena in the regime of 
spontaneous chiral symmetry breaking, two further remarks are important: 

1. at non-zero lattice spacing the twisted mass term nxilb'TzX induces the twisted conden- 
sate ([x^75T3x]r), while the untwisted mass terms x^\U]x of sq. (jT^ determines the 
untwisted condensate ([xx]r)- 

2. the local plaquette field 

'^(^) ^hX- \tT[U,{x)U,{x + afi)Ul{x + aP)Ul{x)] (15) 

admits on the basis of lattice symmetries an operator expansion of the form 

(/)(x) = 6ol + hg a^F ■ F + 63 a^[xx]sub + &4 o'^/iix^Ts^axlsub + O(a^) , (16) 

with [...]sub denoting a subtracted, multiplicatively renormalizable, operator and F the 
continuum gauge field strength tensor. The plaquette expectation value P{r, aniQ, a^) 
can be correspondingly written in the form 

P{r,aniQ,afi) = hQ + higa^{F ■ F)^r,amo,ap) + ^3a^([XX]sub)(r,amo,a/.) 

+ &4aV([X^75^3X]sub)(r,amo,a/.) + O(a^) . (17) 

The important point about the representation p7|) is that it shows that P is actually 
sensitive to the value of the subtracted condensates ([x^75T3x]sub) and ([xx]sub)- 

Before continuing it is useful to pause a moment and discuss the structure of eqs. (fTB|) and (fT7|) 
and nature of the various terms appearing in it. 
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We first notice that the contributions from the identity and the F ■ F operator are 
put together within a square parenthesis in eqs. (pi)|) - |T7j) to remind us that there is 
no unambiguous way to subtract from the latter its power divergent mixing with the 
identity. Ultimately this is due to the fact that, unlike the chiral condensates, the 
vacuum expectation value oi F ■ F is not an order parameter of any symmetry. 

For the reason we have just recalled, it is instead perfectly possible to unambigously 
define, in the massless limit, multiplicatively renormalizable operators [xx]sub and 
[x^75''"3X]sub5 by following the procedure outlined in refs. |2I]. More generally, such 
quark bilinears can be defined as finite operators even at non-vanishing masses, though 
not uniquely. This can be done by setting, for instance 

[xx]sub = XX- a'^Cso (r, am, a/j,) , (18) 

[x«75^3X]sub = xhbTsX - a" V Cp{r, am, a^i) , (19) 

with the dimensionless coefficient functions Cgo and Cp determined at some finite space- 
time volume l^ = Vo by the conditions 

([XX]sub)(r,™o,/.) = 0, V = V^, (20) 

([X^75T-3X]sub)(r,mo,M) = , V = Vq . (21) 

Both the coefficients C50 and Cp admit a finite polynomial expansion in afh and a/i 
(actually in {a^Y fo^ parity reasons). 

In terms of [xx]sub and [x^75T3x]sub, the renormalized scalar density in the physical 
basis, [ipip\R, reads 

[^iIj]r = Z^l{uj)Zp[cosuj[xx]suh + sinu;[xi75r3X]sub] , (22) 

where Zm = Zp/Zgo, Zu = [-s^mCos^cu + sin^o;]^/^ and Zr denotes the renormalization 
constant of xTx in the standard Wilson regularization computed in a mass independent 
renormalization scheme ^. Consistently with the general arguments given above, we re- 
mark that only the leading a~^ divergent subtraction is uniquely fixed by the symmetries 
of the theory (WTI's and spurionic transformations). Consequently these properties 
can be used to make the chiral scalar condensate, ipip, multiplicative renormalizable in 
the massless limit, by defining it, e.g. as the Wilson average over the expectation values 
computed with opposite values of the coefficient of the Wilson term [221 E] • 



^Thc relations between renormalized and subtracted operators in the ^-basis are [xx]r ~ Zs"[xx]suh and 
[xil5T3X\R = Zpixi-fsTaxhuh- 
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After this little digression let us go back and discuss the implications of eq. (jl7|l . If we are 
on the lattice and take the action of eq. (jH)) for values of /i or m much larger than 0(aAQQj-,), 
the condensates ([x^TsTsxlsub) or ([xx]sub) are expected to show no metastability and thus 
the same should be true for the plaquette expectation value. However, if fi is smaller than 
OlaA^QY)) the physical scalar condensate signaling spontaneous chiral symmetry breaking is 
not simply given by ([x^Ts^sxlsub), but has in general also an untwisted component, ([xx]sub)- 
Both components have an impact on the value of the plaquette (see eq. fll7|) ). When m passes 
from positive to negative values the expectation value of the untwisted operator [xx]sub 
should also change sign and, at non-vanishingly small values of fi, eventually become very 
small for almost critical values of rriQ. In this situation, owing to the presence of the chiral 
symmetry breaking term ()13p in the action, the tmQCD sample of gauge configurations 
is expected to include configurations where {[xx]suh)u is positive and configurations where 
([xx]sub)(7 is negative, corresponding to whether m^^^^ is positive or negative, respectively. 
(By (...){/ we mean the fermionic Wick contraction on a fixed gauge background U.) Since 
the coefficient 63 = 63 (r, arriQ, afi) does not vanish at ttiq = "icrit ^, the value of the plaquette 
on the configurations where ([xx]sub)t/ is positive should be different - on the basis of the 
operator expansion (fT^ - from that on the configurations where ([xx]sub)(7 is negative. The 
observed jumps of the plaquette expectation value can hence be regarded as a combined 
effect of spontaneous chiral symmetry breaking and the explicit breaking of this symmetry 
due to the Wilson term in eq. (fT^ . 

4.2 Effective potential model 

The scenario of a jump in the scalar condensate for Wilson fermions on the lattice has 
actually been given already some time ago by Sharpe and Singleton [22 • As it has been 
shown in that work, the phase structure of lattice QCD for /i = with Wilson-type quarks 
can be understood in the low energy chiral theory of pseudo-Goldstone bosons if the influence 
of leading lattice artifacts of 0{a) and O(a^) is taken into account. 

There are two alternatives: either there exists an Aoki phase |22] or there is a flrst order 
phase transition between the phases with positive and negative quark mass and the Aoki 
phase does not exist. 

The relevant part of the effective potential is written in [22] as 

V^ = -ciA + C2A^ . (23) 



■^Using the spurionic invariances of the action ^, it is possible to show that 63 is odd under (r -^ 
— ?') X (toq -^ —Too), or equivalently, since md—r) — —Tnc{r), under (r -^ — r) x (m -^ — m). We expect 
hence a contribution to 63 odd in r and even in m. 



Here A denotes the flavour singlet component of the SU(2) matrix valued field S in the low 
energy effective chiral Lagrangian: 
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j: = A + iY^ BrTr . (24) 



Because of the relation 1 = A^ + X]r=i BrBr the variable A is constrained to lie between -1 
and +1 inclusive. In the vicinity of the critical quark mass the constant C2 = C^(a^) and the 
other parameter Ci is proportional to the bare quark mass (in our notations Ci oc m). 

In order to find the ground state ( "vacuum" ) the effective potential has to be minimized. 
Without repeating the details of the discussion in ref. [22] let us just summarize the result. 

In case of positive C2 there exists an Aoki phase in the region of bare quark masses 
defined by — 2c2 < Ci < 2c2. At the boundaries ci = ±2c2 all three pion masses vanish. 
Inside the Aoki phase the charged pions are massless because they are the Goldstone bosons 
of spontaneous flavour symmetry breaking but the neutral pion is massive. Outside the 
Aoki phase (|ci| > 2c2) the flavour symmetry is preserved by the ground state and the three 
degenerate pions are massive (see figure IHl). 

The other alternative is that C2 is negative. In this case the flavour symmetry is preserved 
everywhere but there exists a minimal pion mass because the pion mass is given by 

ml = f^'ilc^l + 2\C2\) . (25) 

At ci = the vacuum expectation value jumps from H = A = +1 to T, = A = —1. Since 
the jump of this "order parameter" happens at non-zero pion mass (i.e. finite correlation 
length) the thermodynamical description of the behaviour near ci = corresponds to a first 
order phase transition. 

An interesting intermediate situation is defined by C2 = 0. In this case the vacuum 
expectation value jumps between S = A = +1 and H = A = —1 at a. single first order 
phase transition point. This limiting case is the ideal situation, when the phase structure 
in the Sharpe-Singleton model is identical to the expected one in the continuum. It can be 
characterized either by saying that the Aoki phase has zero extension or that the minimal 
pion mass is zero (see figure Ej). Of course, this behaviour is valid only up to 0{a^) effects, 
neglecting higher orders in the chiral expansion. 

4.3 Scenarios 

Our numerical results reveal that we clearly observe metastabilities in various quantities. 
Thus our conclusion is that at least for vanishing twisted mass parameter, i.e. for the standard 
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Wilson lattice theory, there is a first order phase transition. For non- vanishing values of /x 
we can have two scenarios. 

The first is that the first order phase transition persists for /i 7^ but sufficiently small in 
absolute value. For large fi the theory approaches the quenched limit with a constant quark 
determinant and therefore it is plausible that no phase transition is expected. This scenario 
suggests that the first order phase transition line in the (mo, A*) plane has an end point: the 
two phases with positive and negative quark masses are analytically connected (see figure 
ini). The situation is in this sense analogous to the phase structure of the SU(2) fundamental 
Higgs model (see chapter 6 of J23] and references therein). 

The second scenario is that for any non- vanishing value of fi the first order phase transi- 
tion disappears. In this scenario, when varying mo, one passes at some small distance from 
the first order phase transition at // = and just feels this close- by phase transition. 

We can at present not differentiate between these two scenarios. From the numerical 
side we would need to know better the /i and /? dependence of the metastability phenomena. 
From the analytical side an analysis a la Sharpe and Singleton including the twisted mass 
parameter /i is helpful.^ 

The first order phase transition between the phases with positive and negative quark 
masses observed in the previous section is consistent with the no-Aoki-phase alternative 
(c2 < 0) of Sharpe and Singleton. 

Our exclusion of the Aoki phase is in agreement with the results of a recent paper [23] 
which suggests that in case of the unimproved Wilson action the Aoki phase is restricted 
to the region of strong gauge couplings (/5 < 4.6). Note that in an early paper on QCD 
thermodynamics with Wilson quarks jSHj a first order "bulk" phase transition has also been 
observed at /5 = 4.8 which is consistent both with ref. J2S1 and with our observations. For 
further numerical work on the Aoki phase, see refs. fI7\ 

The rather strong metastability of the two phases with positive and negative quark mass 
can be understood on the basis of the properties of the eigenvalue spectrum of the (non- 
hermitean) Wilson-fermion matrix in the twisted mass basis corresponding to eq. (jH)). For 
zero twisted mass (/i = 0) at small positive quark masses there is a "cloud" of eigenvalues 
close to the origin near the real axis. (For a numerical study see section 4 of J7].) In order 
to reach negative quark masses this "cloud" has to cross near the origin to the other side 
with negative real parts. This eigenvalue cloud crossing is strongly suppressed by the zero of 
the determinant. This, we believe, is the reason at the microscopical level for the observed 
strong metastability. For non-zero twisted mass there is a strip of width 2|/i| around the 
real axis where there are no eigenvalues. If this strip is wide enough the eigenvalues are 



^Wc thank Gemot Miinster for discussions on this and for communicating us his results before pubhcation. 
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Figure 6: The alternatives of the phase structure in the {rriQ, fi)- 
plane: Upper part: Aoki phase at /i = if C2 > 0, middle part: 
first order phase transition point if C2 = 0, lower part: first order 
phase transition hue if C2 < 0. In the latter case the two phases are 
connected with each other as it is shown by the curve with arrows 
at both ends. 



21 



sufficiently far away from the origin and the first order phase transition disappears. 

As it was aheady emphasized in |2S] , the sign of the coefficient C2 in the low energy pion 
effective potential is not universal, it depends on the way the action is discretized. Therefore 
a clever choice of the lattice action may weaken the ffist order phase transition and, for 
instance, decrease the minimal pion mass at it. Previous results of the JLQCD Collaboration 
fK)\ support the conjecture that changing the gauge action alone has an important effect. 
If, indeed, one could find some parameter in the lattice gauge action which at some value 
would change the sign of C2 an appealing possibility would be to tune the lattice action to 
this value. The features of a discretization with C2 = seem to be quite favourable from 
the point of view of light quark simulations when, up to C(a^), there would be just a single 
point in the {rriQ^fi) plane with vanishing pion mass - an ideal situation corresponding to 
the expected phase structure in the continuum. 

5 Conclusion 

In this paper we have explored Wilson twisted mass fermions restricting ourselves to simula- 
tions at only one value of /? = 5.2. By fixing the twisted mass parameter /i and changing the 
untwisted Wilson quark mass mo, or equivalently the hopping parameter k, we encountered 
strong metastabilities in the plaquette expectation value, visible both in thermal cycles as 
well as in long-living metastable states. At the same time, the pion mass does not vanish 
but has a minimum at a rather large value. The PCAC quark mass m^^'^'^ in the different 
metastable branches is positive for the branch with low plaquette expectation value and it 
is negative for the branch with high plaquette expectation value. 

The detection of these metastabilities became possible by employing a twisted mass 
term. Only a non-vanishing value of /i allowed us to cross the critical quark mass. We 
showed that for lattice theories that break chiral symmetry explicitly the jump of the scalar 
condensate, when changing the sign of the quark mass, induces a jump of the plaquette 
expectation value with associated signs of metastability. For /i = these metastabilities 
find a natural interpretation in the effective potential model of Sharpe and Singleton, arising 
from spontaneous symmetry breaking and using a low energy effective Lagrangian which also 
describes lattice artifacts. The agreement with the Sharpe-Singleton model is remarkable 
because in the continuum limit in this model the phase structure of lattice QCD with Wilson 
quarks approaches fast - at a rate O(a^) - the expected phase structure of QCD near 
zero quark mass. This is an important property which has to be required from any lattice 
regularization of QCD. 

It should be clear that our work can only represent a first step in a detailed understanding 



22 



of the QCD phase diagram at zero temperature near vanishing quark masses. Clearly, 
substantially more work has to be done to resolve this phase structure and its behaviour 
in the continuum limit. For instance, at present for /_i 7^ we are unable to differentiate 
between a scenario where the first order phase transition persists and another one where at 
yU 7^ only a remnant of the phase transition at /i = is seen. In this respect an analysis 
like in ref. [SB] for /i 7^ is very helpful |28j . 

Among the many open questions there are: How fast does the gap vanish when the 
continuum limit at higher values of (3 is approached? How are the signs of metastability 
related to the ones observed using the Wilson plaquette action and clover-improved Wilson 
fermions? How precisely do the eigenvalues re-arrange when the critical quark mass is 
crossed? Do different gauge actions change the couplings of the effective potential and may 
hence lead to avoid the phenomena of metastability and reproduce the ideal phase structure 
at vanishing quark mass already for non-zero lattice spacing? 

The most important question is, of course, how phenomenology can be done given the 
metastability phenomenon seen in our present results; i.e. what is the lowest value of the 
quark mass that can be reached before one enters the regime of metastabilities and how does 
this change with decreasing value of the lattice spacing. 
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A Appendix 

A.l Even-odd preconditioning for the HMC algorithm 

Let us start with the Dirac operator in the hopping parameter representation in the twisted 
basis written as 

5[x,X,t/]=$^x(a:)M.,x(l/) (26) 

xy 

where the matrix M can be easily read from eq. ((Tj). Using M one can define the hermitian 
operator 

where the submatrices Q± can be factorised as follows: 



eo 



'75M± 0\ (l {Mtr'M,, 



(2^ 



and we have defined jl = 2k^. Note that (M^)~^ can be easily computed to 

Using det((5) = det((5+) det((5_) one can now derive the following relation (an equation 
apart from an irrelevant factor): 

det((5±) oc det((0±) 

(29) 
Q± ■■= l^{Mi - M„e(M±)-iM,,) , 

where Q± is only defined on the odd sites of the lattice. In the HMC algorithm the deter- 
minant is stochastically estimated using pseudo-fermion fields (po'- 

det(g+g_) = j Z}[0o, 0l] exp(-5fc) , 

5,:=0l(g+g_)-Vo, 

where the fields 0o are defined only on the odd sites of the lattice. In order to compute the 
force corresponding to the effective action Si, we need the variation of S^ with respect to the 
gauge fields (using 5{A'^) = —A^^5AA~^): 

(30) 
= -[Xl6Q+Yo + Y^6Q^Xo] 
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with Xo and Yq defined on the odd sites as 

X, = {Q+Q^)-'<p,, Yo = Q~'<p, = g„X„ , (31) 

where (5± = Q^f has been used. The variation of Q± reads 

SQ± = 75 (-5M„,(M±)-iM,„ - MUM^,)-'6M,,) , (32) 

and one finds 

SSb = -{X^6Q+Y + Y^SQ^X) = -{X^6Q+Y + (X^SQ+YY) (33) 

where X, Y is now defined over the full lattice as 

^ ^ /-(M-)-iMeoX,\ ^ ^ ^ /-(M+)"iM,,n\ _ ^g^^ 

In addition, 6Q+ = SQ-, Mj^ = ■y^Moej^ and M^g = 75Meo75 has been used. Since the 
bosonic part is quadratic in the 0o fields, the 0o are generated at the beginning of each 
molecular dynamics trajectory with 

(po = Q+R, (35) 

where i? is a random spinor field taken from a Gaussian distribution with norm one. 

A. 1.1 Hasenbusch trick 

The trick first presented in ^U] is based on the observation that writing 

det[Q+g_] = det[W+W.] ■ det[iQ+Q.)/{W+W.)] (36) 

is advantageous for the HMC, if the condition number of W+W_ and of {Q^QJ)/{W^WJ) is 
significantly reduced compared to the condition number of only {Qj^QJ). In order to achieve 
this we define 

Q± = l^Dw ± i/i , 

(37) 

W± = 75-Diy ± «/i2 • 

With /i2 = /i + A/i it follows immediately that the condition number of W+W- is lower 
than the one of Q+Q^ if for Amin and A^ax the lowest and the largest eigenvalue of Q+Q^, 
respectively, |Ajnin| <ti fi^ '^ l^maxl holds: the condition number of W^W_ is lAmaxl/yui while 
the one of (W+W-)~^{Q-^^Q^y contrariwise is /il/IAmml- We can take fi which is a lower 
bound for |Amin| to write down the following estimates for the condition numbers k: 

I \ I r/2 

, l^maxl ; ^ P2 

f«W+W. - ~2 ' fc(Q+Q-)/(W+W_) < -T^ , 
P2 H- 
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iP^i-~^')iQ+Q^r'. 



which leads to an optimal choice for p,2 = v^l-^maxl ■ At^- As has been shown in ^T] also the 
force contribution coming from {Q+Q_)/(W+W-) is reduced. This is true also for tmQCD 
and can bee seen in the following way: noticing that 

Q+Q_ = Q'^ + jl^ and 

it follows that 

W+W_{Q+Q^)-' = 1 

Since the corresponding effective action reads 

one can see that one gets an explicit factor {fi\ — /i^) <^ 1 multiplying the force contribution 
compared to the original effective action which will reduce the force and therefore lead to a 
smoother evolution of the algorithm. 

Let us remark that the procedure explained above can be immediately applied to the 
even-odd preconditioned system. Furthermore the trick can be iterated to two or even more 
additional operators. 

5500 



- 5000 



(38) 

(39) 
(40) 



Pa 




- 4500 



C 



- 4000 



- 3500 



3000 



Figure 7: Acceptance rate Pa and cost C in units of CG iterations versus /i2 = ^^^2 at fixed 
HMC stepsize and trajectory length. The dashed line represents the cost required to obtain 
about 90% acceptance rate without the additional operator. The parameters are: 8^ lattice, 
/3 = 5.2, /€ = 0.17, /i = 0.01. 



In fig. [7| the cost C in units of CG iterations and the acceptance rate Pa is plotted versus 
/i2 = 2hi^2 at fixed HMC stepsize and trajectory length. One can see that as expected the 
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acceptance rate increases by introducing an additional operator and reaches a maximum 
around fi = 0.2. Of course also the costs increase when compared to the HMC without 
additional operators. But the costs are still much less than what is needed to reach an 
acceptance rate of about 90% without the additional operator (see the dashed line in fig. E)). 
One can see that the gain in the costs is about a factor of two. 

A. 2 Even-odd preconditioning for the TSMB algorithm 

In this appendix even-odd preconditioning is derived for the TSMB algorithm. The even-odd 
subspace decomposition of the fermion matrix in the twisted basis can be written as 

y -\Moe /il +«75^3/W j 

where indices start by zero = even, the lattice spacing is set to a = 1 and the abbreviation 
jii = rriQ + Ar = (2k)~^ is introduced. The hermitean fermion matrix Q = •y^TiQ^ = Q^ is 
then 

\ -\l->TiMoe 1-,Tilli + T2H J 

Using the notation 

^5 = (75n/ii + i"2/i)~S57-i = (/il - i75r3/i)(/ii + i?y^ (43) 

one can write Q as the following product: 



Q 



i^Tijii + T2H W 1 

75nAil + r2/i j I -^t^Moe 1 



1 \ / 1 -IhMeo 



1-^hMoehMe, 




(44) 



This can be used for preconditioned inversion of Q because the inverse of all the factors 
but the third one is trivial. Of course, the third factor is expected to have smaller condition 
number than Q itself. 

Multi-boson (MB) updating can be set up following [TH] . Since the determinant of the 
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above triangular matrices is equal to 1 we have 

detg = detl ^^-^/^^ + -^^ ^0 

= det (75ri/ii + T2fi) 

e 

■ det ( 75ri/ii + T2fi - -75nMoe(75ri/ii + r2/i)"S5nMeo J (45) 

where dete and deto denote determinants in the even and odd subspaces, respectively. The 
first factor does not depend on the gauge field and therefore it can be omitted. In the second 
factor we have the hermitean matrix defined on odd sites 

Q = 75T1/XI + r2/i - -75TlMoe(75ri/ii + r2/i)"S5TlMeo = 

= 75^iAH + T-2/i - T75nMoe(75ri/ii + r2/i)(/ii + /i^)"S5nMeo = (Q"^ . (46) 

The hermiticity of Q, which can be called hermitean preconditioned fermion matrix, follows 
from 

75riM],75ri = M,„ . (47) 

In MB updating one can start with the identity 

det Q = det g oc det g = f det qA ' ~ -—^— (48) 

o V o / deto_Pi(g^) 

where the Pi is a polynomial approximation satisfying 

P.{x)^\ (49) 

in an interval x G [e, A] covering the spectrum of Q^. (Note that for fi ^ det Q and det Q 
are positive.) 

The rest is the same as usual: one writes the polynomial with the help of the square 
roots of its roots pj, j = 1, 2, . . . as 

P^m cxn(0-p*)(g-p,). (50) 

j 

Then using the identity 

det I /' /° ) = det A,, ■ det {A^o - A,,A-^'A,,) (51) 



A A e 
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one obtains 



det(Q - pj) = det (75X1/11 + r2/i)" 



det 



-^-f^TiMoe 75n/il + 7-2/i - pj 



(52) 



Denoting the projector on the odd subspace by Po we finally obtain the multi-boson repre- 
sentation 



detg')' oc Yl 



1 



det 



(Q - PoP*)iQ - PoPj) 



oc 



jm exp I - J] $](g - P,p*){Q - P,p,)^, 



(53) 
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